module shell_mod
    use iso_fortran_env, only: real64,int64
    implicit none
    !> 单元自由度数
    integer, parameter :: elem_dof_num=24
    !> 节点自由度数
    integer, parameter :: node_dof_num=6
    !> 单元节点数
    integer, parameter :: el_node_num=4
    !> 四个节点的自然坐标(xi,eta)
    real(real64), parameter :: node_natural_coords(4,2) = reshape([ 1.0_real64,-1.0_real64,-1.0_real64, 1.0_real64,&
                                                                    1.0_real64, 1.0_real64,-1.0_real64,-1.0_real64], [4, 2])
    ! 四节点单元的高斯积分点位置 
    real(real64), parameter    :: p_=sqrt(1.0_real64/3.0_real64)
    !> 高斯积分权重数组
    real(real64),parameter     :: gauss_weight(4) = [1.0_real64, 1.0_real64, 1.0_real64, 1.0_real64]
    !> 高斯积分等参xi坐标数组
    real(real64),parameter     :: gauss_xi_coord(4) = [p_,-p_,-p_,p_]
    !> 高斯积分等参eta坐标数组
    real(real64),parameter     :: gauss_eta_coord(4) = [p_,p_,-p_,-p_]
    !> 板的横向剪切系数
    real(real64),parameter     :: k_=5.0_real64/6.0_real64
    !> 四节点壳单元的膜刚度矩阵相关自由度的索引
    integer, parameter :: kme_inds(12)=[1,2,6,7,8,12,13,14,18,19,20,24]
    !> 四节点壳单元的板弯曲和剪切刚度矩阵相关自由度的索引
    integer, parameter :: kbs_inds(12)=[3,4,5,9,10,11,15,16,17,21,22,23]
contains
! ---------------------------------------------------------
! 计算平面壳的单元刚度矩阵ke=kme+(kbe+kse) (单元局部坐标系下)
! ---------------------------------------------------------
    !> 计算第i个节点的第dofth自由度的全局自由度
    function get_dof_ind(i,dofth) result(dof_ind)
        integer,intent(in):: i,dofth
        integer::dof_ind
        dof_ind = i*node_dof_num-(node_dof_num-dofth)
    end function get_dof_ind

    !> 计算节点的自由度索引数组(3)
    function get_ndof_inds(i) result(retval)
        integer, intent(in) :: i
        integer :: retval(6)
        retval(1) = 6*i-5
        retval(2) = 6*i-4
        retval(3) = 6*i-3
        retval(4) = 6*i-2
        retval(5) = 6*i-1
        retval(6) = 6*i
    end function get_ndof_inds

    !> 计算单元的全局自由度索引数组
    subroutine get_eldof_inds(el_conn,el_inds)
        implicit none
        integer, intent(in) :: el_conn(4)
        integer, intent(out) ::  el_inds(24)
        integer:: i,j
        el_inds(:)=0
        do i=1,4
            j=el_conn(i)
            el_inds(6*i-5)=6*j-5
            el_inds(6*i-4)=6*j-4
            el_inds(6*i-3)=6*j-3
            el_inds(6*i-2)=6*j-2
            el_inds(6*i-1)=6*j-1
            el_inds(6*i)=6*j
        end do
    end subroutine get_eldof_inds
    
    !> 获取壳单元在全局坐标系下的单元集中质量矩阵(24,24)
    subroutine get_gme(rho,t,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,gme)
        ! 密度,厚度,全局节点坐标
        real(real64),intent(in):: rho,t,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4
        real(real64),intent(out):: gme(24,24)
        integer:: i
        real(real64):: area,a,b,c,d
        gme(:,:)=0.0_real64
        do i=1,24
            gme(i,i)=1.0_real64
        end do 

        do i=1,4
            gme(6*i-2,6*i-2)=gme(6*i-2,6*i-2)*((t*t)/12.0_real64)
            gme(6*i-1,6*i-1)=gme(6*i-1,6*i-1)*((t*t)/12.0_real64)
            gme(6*i,6*i)=gme(6*i,6*i)*((t*t)/12.0_real64)
        end do

        ! 计算单元面积
        a=sqrt((gx2-gx1)*(gx2-gx1)+(gy2-gy1)*(gy2-gy1)+(gz2-gz1)*(gz2-gz1))
        b=sqrt((gx3-gx2)*(gx3-gx2)+(gy3-gy2)*(gy3-gy2)+(gz3-gz2)*(gz3-gz2))
        c=sqrt((gx4-gx3)*(gx4-gx3)+(gy4-gy3)*(gy4-gy3)+(gz4-gz3)*(gz4-gz3))
        d=sqrt((gx1-gx4)*(gx1-gx4)+(gy1-gy4)*(gy1-gy4)+(gz1-gz4)*(gz1-gz4))
        area=0.25*sqrt((-a+b+c+d)*(a-b+c+d)*(a+b-c+d)*(a+b+c-d))

        gme=gme*(rho*area*t)*0.25_real64
    end subroutine get_gme
    
    !> 输入材料属性和单元全局节点坐标(3,4), 获取壳单元在全局坐标系下的单元刚度矩阵gke(24,24)
    subroutine get_gke(e,nu,t,shell_coords,gke)
        implicit none
        real(real64), intent(in) :: e,nu,t,shell_coords(3,4)
        real(real64), intent(out) :: gke(24,24)
        ! internal vars
        real(real64):: gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4
        real(real64):: kme(12,12),kbe(12,12),kse(12,12),tmatrix(24,24)
        real(real64):: sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4 ! 壳单元局部坐标系下的节点坐标
        gke(:,:)=0.0_real64
        kme(:,:)=0.0_real64
        kbe(:,:)=0.0_real64
        kse(:,:)=0.0_real64
        tmatrix(:,:)=0.0_real64
        
        gx1=shell_coords(1,1)
        gx2=shell_coords(1,2)
        gx3=shell_coords(1,3)
        gx4=shell_coords(1,4)
        gy1=shell_coords(2,1)
        gy2=shell_coords(2,2)
        gy3=shell_coords(2,3)
        gy4=shell_coords(2,4)
        gz1=shell_coords(3,1)
        gz2=shell_coords(3,2)
        gz3=shell_coords(3,3)
        gz4=shell_coords(3,4)

        ! 计算节点在壳单元局部坐标系下的坐标
        call get_shell_sys_coords(gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)
        
        ! 计算膜刚度矩阵kme
        call get_membrance_ke(e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,kme)
        gke(kme_inds,kme_inds)=gke(kme_inds,kme_inds)+kme  ! 将kme加到gke中

        ! 计算板弯曲刚度矩阵kbe
        call get_bend_ke(e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,kbe)
        gke(kbs_inds,kbs_inds)=gke(kbs_inds,kbs_inds)+kbe  ! 将kbe加到gke中
        
        ! 计算剪切刚度矩阵kse
        call get_shear_ke(e,nu,t,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,kse)
        gke(kbs_inds,kbs_inds)=gke(kbs_inds,kbs_inds)+kse  ! 将kse加到gke中

        call get_tmatrix(gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,tmatrix)
        gke=matmul(transpose(tmatrix),matmul(gke,tmatrix))  ! 将tmatrix作用到gke上
    end subroutine get_gke

    !> 单元刚度矩阵kbe对应自由度(12,1): {w1_0,theta1_x,theta1_y,...,w4_0,theta4_x,theta4_y}^T
    subroutine get_bend_ke(e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,kbe)
        implicit none
        real(real64), intent(in) :: e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4
        real(real64), intent(out) :: kbe(12,12)
        ! internal vars
        integer::i
        real(real64)::ixi,ieta,d(3,3),bb(3,12),detj
        
        ! 初始化
        kbe(:,:)=0.0_real64
        ! 高斯积分
        d=d_bend(e,t,nu)
        do i=1,4
            ixi=gauss_xi_coord(i)
            ieta=gauss_eta_coord(i)
            call get_bend_bmatrix(ixi,ieta,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,bb)
            detj=det_jacobian(ixi,ieta,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)
            kbe=kbe+gauss_weight(i)*t*detj*matmul(transpose(bb),matmul(d,bb))
        end do
    end subroutine get_bend_ke

    !> 单元刚度矩阵对应的自由度(12,1): {w1,theta1_x,theta1_y,...,w4,theta4_x,theta4_y}^T
    subroutine get_shear_ke(e,nu,t,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4&
                            ,gz1,gz2,gz3,gz4,kse)
        ! > @param : 使用混合插值方法构造板的剪切刚度矩阵kse
        implicit none
        real(real64), intent(in) :: e,nu,t,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4
        real(real64), intent(out) :: kse(12,12)
        ! internal vars
        real(real64)::ixi,ieta,d(2,2),bs(2,12),global_p(3),local_p(3),shell_sys_origin(3),&
                shell_sys_x(3),shell_sys_y(3),shell_sys_z(3)
        real(real64)::sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4 ! 壳单元局部坐标系下的节点坐标
        integer::i

        ! 初始化
        d(:,:)   = 0.0_real64
        kse(:,:) = 0.0_real64
        bs(:,:)  = 0.0_real64
        local_p(:)=0.0_real64
        shell_sys_origin(:)=0.0_real64
        shell_sys_x(:)=0.0_real64
        shell_sys_y(:)=0.0_real64
        shell_sys_z(:)=0.0_real64

        ! 计算节点在壳单元局部坐标系下的坐标值
        shell_sys_origin=0.25_real64*[gx1+gx2+gx3+gx4, gy1+gy2+gy3+gy4, gz1+gz2+gz3+gz4]
        shell_sys_x=0.5_real64*[gx1+gx4, gy1+gy4, gz1+gz4]-shell_sys_origin
        shell_sys_y=0.5_real64*[gx1+gx2, gy1+gy2, gz1+gz4]-shell_sys_origin
        shell_sys_z=cross_product(shell_sys_x,shell_sys_y)
        shell_sys_x=shell_sys_x/norm2(shell_sys_x)
        shell_sys_y=shell_sys_y/norm2(shell_sys_y)
        shell_sys_z=shell_sys_z/norm2(shell_sys_z)
        
        global_p=[gx1,gy1,gz1]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx1=local_p(1)
        sy1=local_p(2)
        
        global_p=[gx2,gy2,gz2]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx2=local_p(1)
        sy2=local_p(2)
        
        global_p=[gx3,gy3,gz3]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx3=local_p(1)
        sy3=local_p(2)

        global_p=[gx4,gy4,gz4]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx4=local_p(1)
        sy4=local_p(2)
        
        ! 单元的剪切应力-应变矩阵
        d=d_shear(e,nu,t)
        ! 四点高斯积分
        do i=1,4
            ixi=gauss_xi_coord(i)
            ieta=gauss_eta_coord(i)
            call get_shear_bmatrix(ixi,ieta,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,&
                                    sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,bs)
            kse=kse+gauss_weight(i)*t*det_jacobian(ixi,ieta,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)*matmul(transpose(bs),matmul(d,bs))
        end do
    end subroutine get_shear_ke
    
    !> 单元刚度矩阵kme对应自由度(12,1): {u1_0,v1_0,theta1_z,...,u4_0,v4_0,theta4_z}^T
    subroutine get_membrance_ke(e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,kme)
        !> 计算膜单元的单元刚度矩阵(12,12),与节点的平动自由度u,v以及节点附加旋转角度theta_z有关
        ! 后续平动自由度部分的膜单元刚度矩阵需要采用EAS方法修正
        ! ke=[kuu,kut;
        !     ktu,ktt]
        use matrix_utils_mod , only: inv=>invert_matrix_d ! 重命名,将inv指向invert_matrix
        implicit none
        real(real64), intent(in) :: e,nu,t,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4
        real(real64), intent(out) :: kme(12,12)
        real(real64)::bmu(3,8),bmt(3,4),d_matrix(3,3)
        !! 将膜单元刚度的子矩阵, u:和u,v自由度有相关, t:和theta_z自由度有关
        real(real64)::kuu(8,8),ktt(4,4),kut(8,4),ktu(4,8) 
        !! EAS修正时, 和内部自由度的应变矩阵b_a,刚度子矩阵; a:和内部自由度有关, u:和u,v自由度有关
        real(real64)::b_a(3,5),kaa(5,5),kua(8,5),kau(5,8),kaa_inv(5,5)
        integer::i,j
        real(real64)::t1
        ! 初始化
        i=0
        j=0
        t1=0.0_real64
        bmu(:,:)=0.0_real64
        bmt(:,:)=0.0_real64
        kuu(:,:)=0.0_real64
        ktt(:,:)=0.0_real64
        kut(:,:)=0.0_real64
        ktu(:,:)=0.0_real64
        kme(:,:)=0.0_real64
        b_a(:,:)=0.0_real64
        kaa(:,:)=0.0_real64
        kua(:,:)=0.0_real64
        kau(:,:)=0.0_real64
        d_matrix(:,:)=0.0_real64
        kaa_inv(:,:)=0.0_real64
        ! 第一步,先用Bm矩阵求出具有附加旋转自由度的膜单元刚度矩阵,kuu,ktt..etc,此时未EAS修正.
        d_matrix=dm(e,nu) 
        do i=1,4
            t1=det_jacobian(gauss_xi_coord(i),gauss_eta_coord(i),sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)
            call get_elstrain_matrix(gauss_xi_coord(i),gauss_eta_coord(i),sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,bmu,bmt)
            kuu=kuu+gauss_weight(i)*t*t1*(matmul(transpose(bmu),matmul(d_matrix,bmu)))
            ktt=ktt+gauss_weight(i)*t*t1*(matmul(transpose(bmt),matmul(d_matrix,bmt)))
            kut=kut+gauss_weight(i)*t*t1*(matmul(transpose(bmu),matmul(d_matrix,bmt)))
            ktu=ktu+gauss_weight(i)*t*t1*(matmul(transpose(bmt),matmul(d_matrix,bmu)))
        end do
        
        ! 第二步,使用EAS修正平动自由度对应的单元刚度子矩阵kuu, 从而解决膜内剪切自锁问题
        ! 2.1 计算内部自由度的应变矩阵b_a,刚度子矩阵kaa,kua,kau
        do i=1,4
            t1=det_jacobian(gauss_xi_coord(i),gauss_eta_coord(i),sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)
            ! 计算内部自由度的应变矩阵b_a
            call get_elstrain_matrix_eas(gauss_xi_coord(i),gauss_eta_coord(i),sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,b_a)
            kaa=kaa+gauss_weight(i)*t*t1*(matmul(transpose(b_a),matmul(d_matrix,b_a)))
            kua=kua+gauss_weight(i)*t*t1*(matmul(transpose(bmu),matmul(d_matrix,b_a)))
            kau=kau+gauss_weight(i)*t*t1*(matmul(transpose(b_a),matmul(d_matrix,bmu)))
        end do
        
        ! 2.2 使用静力凝聚技术,计算EAS修正后的刚度矩阵kuu
        call inv(kaa,kaa_inv,j)
        if ( j /= 0 )then
            print*, 'Error in invert_matrix of subroutine : get_membrance_ke'
            stop
        end if
        kuu=kuu-matmul(kua,matmul(kaa_inv,kau)) !修正后的刚度矩阵kuu

        ! 第三步,将刚度矩阵kuu,ktt,kut,ktu组装成单元刚度矩阵kme
        kme([1,2,4,5,7,8,10,11],[1,2,4,5,7,8,10,11])=kuu
        kme([3,6,9,12],[3,6,9,12])=ktt
        kme([1,2,4,5,7,8,10,11],[3,6,9,12])=kut
        kme([3,6,9,12],[1,2,4,5,7,8,10,11])=ktu
    end subroutine get_membrance_ke

! ---------------------------------------------------------
! 辅助子程序/函数,包括应变矩阵,应力应变矩阵,形函数.etc
! ---------------------------------------------------------
    !> 计算膜单元的应变矩阵,(3,12),在局部坐标xi,eta下的矩阵值
    subroutine get_elstrain_matrix(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4,bmu,bmt)
        ! bm矩阵满足{epsilon_xx, epsilon_yy, gamma_xy}^T=bm*{u1_0,v1_0,theta1_z...,u4_0,v4_0,theta4_z}^T
        ! bm=[bmu, bmt] , bmu(3,8), bmt(3,4)
        ! bmu是和节点自由度u，v有关的应变矩阵
        ! bmt是和节点附加旋转角度theta_z有关的应变矩阵
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4,y1,y2,y3,y4
        real(real64),intent(out) ::  bmu(3,8),bmt(3,4)
        real(real64) ::  bm(3,12)
        integer :: i,ipiv(2),info
        real(real64) :: jacb(2,2),b(2)
        b=0.0_real64
        bm=0.0_real64
        bmu=0.0_real64
        bmt=0.0_real64
        jacb=0.0_real64
        ! 遍历四个节点,计算形函数对x,y的偏导数,组装应变矩阵
        do i=1,4
            jacb=f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
            ! 调用lapack的sgesv()求解线性方程组,求解出Ni_x,Ni_y
            b(1)=tshape_func_pxi(i,xi, eta)
            b(2)=tshape_func_peta(i,xi,eta)
            call dgesv(2,1,jacb,2,ipiv,b,2,info)
            if (info /= 0) then
                print*, 'Error in sgesv'
                stop
            end if
            bm(1,3*i-2)=b(1)
            bm(3,3*i-1)=b(1)
            bm(2,3*i-1)=b(2)
            bm(3,3*i-2)=b(2)

            ! 调用lapack的sgesv()求解线性方程组,求解出Ni^{u}_x,Ni^{u}_y
            jacb=f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
            b(1)=u_theta_shape_func_pxi(i,xi, eta,y1,y2,y3,y4)
            b(2)=u_theta_shape_func_peta(i,xi,eta,y1,y2,y3,y4)
            call dgesv(2,1,jacb,2,ipiv,b,2,info)
            if (info /= 0) then
                print*, 'Error in dgesv'
                stop
            end if
            bm(1,3*i)=b(1)
            bm(3,3*i)=bm(3,3*i)+b(2)

            ! 调用lapack的sgesv()求解线性方程组,求解出Ni^{v}_x,Ni^{v}_y
            jacb=f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
            b(1)=v_theta_shape_func_pxi(i,xi, eta,y1,y2,y3,y4)
            b(2)=v_theta_shape_func_peta(i,xi,eta,y1,y2,y3,y4)
            call sgesv(2,1,jacb,2,ipiv,b,2,info)
            if (info /= 0) then
                print*, 'Error in sgesv'
                stop
            end if
            bm(2,3*i)=b(2)
            bm(3,3*i)=bm(3,3*i)+b(1)
        end do
        ! 从bm(3x12)中提取bmu,bmt
        bmu=bm(:,[1,2,4,5,7,8,10,11])
        bmt=bm(:,[3,6,9,12])
    end subroutine get_elstrain_matrix

    !> 计算EAS修正的膜单元增强应变矩阵
    subroutine get_elstrain_matrix_eas(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4,b_alpha)
        ! 计算EAS修正的膜单元增强应变矩阵
        use matrix_utils_mod, only:det2rank,invert_matrix_d
        implicit none
        integer::info
        real(real64), intent(in)  :: xi,eta,x1,x2,x3,x4,y1,y2,y3,y4
        real(real64), intent(out) ::  b_alpha(3,5)
        real(real64) ::  e_alpha(3,5),f0(3,3),f0_inv(3,3),j0(2,2),detj0,detj
        ! 初始化
        b_alpha(:,:)=0.0_real64
        e_alpha(:,:)=0.0_real64
        f0(:,:)=0.0_real64
        f0_inv(:,:)=0.0_real64
        j0(:,:)=0.0_real64
        ! 计算单元中心处的雅可比矩阵的行列式
        j0(1,1)=0.25_real64*(x2+x3-x1-x4)
        j0(2,1)=0.25_real64*(x3+x4-x1-x2)
        j0(1,2)=0.25_real64*(y2+y3-y1-y4)
        j0(2,2)=0.25_real64*(y3+y4-y1-y2)
        detj0=det2rank(j0)
        
        ! 计算(xi,eta)处的雅可比矩阵的行列式
        detj=det_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
        
        ! 组装增强应变的转换矩阵f0
        f0(1,1)=j0(1,1)*j0(1,1)
        f0(2,1)=j0(1,2)*j0(2,1)
        f0(3,1)=j0(1,1)*j0(2,1)
        f0(1,2)=j0(2,1)*j0(1,2)
        f0(2,2)=j0(2,2)*j0(2,2)
        f0(3,2)=j0(1,2)*j0(2,2)
        f0(1,3)=2.0_real64*j0(1,1)*j0(1,2)
        f0(2,3)=2.0_real64*j0(2,1)*j0(2,2)
        f0(3,3)=j0(1,1)*j0(2,2)+j0(1,2)*j0(2,1)
        call invert_matrix_d(f0,f0_inv,info) ! ! 计算f0的逆矩阵
        
        ! 计算自然坐标系下的增强应变场
        e_alpha(:,:)=0.0_real64
        e_alpha(1,1)=xi
        e_alpha(2,2)=eta
        e_alpha(3,3)=xi
        e_alpha(3,4)=eta
        e_alpha(1,5)=xi*eta
        e_alpha(2,5)=-1.0_real64*xi*eta
        e_alpha(3,5)=xi*xi-eta*eta

        ! 计算EAS修正的增强应变场的应变矩阵
        b_alpha=(detj0/detj)*(matmul(transpose(f0_inv),e_alpha))
    end subroutine get_elstrain_matrix_eas

    !> 计算板单元弯曲应变-位移矩阵
    subroutine get_bend_bmatrix(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4,bb)
        ! 公式来自论文公式2.70, p27
        implicit none
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4,y1,y2,y3,y4
        real(real64),intent(out) ::  bb(3,12)
        ! internal vars
        integer :: i,ipiv(2),info
        real(real64):: jacb(2,2),b(2)
        ! 初始化
        bb(:,:)=0.0_real64
        b(:)=0.0_real64
        jacb(:,:)=0.0_real64
        do i=1,4
            jacb=f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
            ! 调用lapack的sgesv()求解线性方程组,求解出Ni_x,Ni_y
            b(1)=tshape_func_pxi(i,xi, eta)
            b(2)=tshape_func_peta(i,xi,eta)
            call dgesv(2,1,jacb,2,ipiv,b,2,info)
            if (info /= 0) then
                print*, 'Error in sgesv, subroutine : get_bend_bmatrix'
                stop
            end if
            bb(1,3*i)=-1.0_real64*b(1)
            bb(2,3*i-1)=b(2)
            bb(3,3*i-1)=b(1)
            bb(3,3*i)=-1.0_real64*b(2)
        end do
    end subroutine get_bend_bmatrix

    !> 计算板单元的剪切-位移矩阵
    subroutine get_shear_bmatrix(xi,eta,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,&
                                sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4,bs)
        !> 公式来自论文公式2.81, p29
        !> 输入的gx1,gy1...,gx4,gy4是单元全局坐标系的坐标值
        !> sx1...是壳单元局部坐标系下的节点坐标值,用来计算雅可比矩阵
        !> 真正带入Bs公式的时候需要的是节点在板单元坐标系下的坐标值x1,y1等
        !> 全局坐标
        implicit none
        real(real64), intent(in) :: xi,eta,gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,&
                            gz1,gz2,gz3,gz4,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4
        real(real64),intent(out) ::  bs(2,12)
        ! internal vars
        real(real64):: sa,ca,sb,cb ! 分别是 sin(alpha), cos(alpha), sin(beta), cos(beta)
        real(real64):: lbx(3),lby(3),lbz(3),a(3),b(3),c(3),d(3),bd_vec(3),ca_vec(3),&
                v34(3),v32(3),g1(3),g2(3),g3(3),t1(2,2),t2(2,12)
        ! 局部板坐标系下的节点坐标
        real(real64)::x1b,x2b,x3b,x4b,y1b,y2b,y3b,y4b
        ! 公式中的Ax,Bx...等
        real(real64)::ax,ay,bx,by,cx,cy
        ! 初始化
        bs(:,:)=0.0_real64
        lbx(:)=0.0_real64
        lby(:)=0.0_real64
        lbz(:)=0.0_real64
        g1(:)=0.0_real64
        g2(:)=0.0_real64
        g3(:)=0.0_real64
        t1(:,:)=0.0_real64
        t2(:,:)=0.0_real64
        ! 计算a,b,c,d的坐标
        a = 0.5_real64*[gx1+gx2, gy1+gy2, gz1+gz2]
        b = 0.5_real64*[gx2+gx3, gy2+gy3, gz2+gz3]
        c = 0.5_real64*[gx3+gx4, gy3+gy4, gz3+gz4]
        d = 0.5_real64*[gx4+gx1, gy4+gy1, gz4+gz1]
        bd_vec=[b(1)-d(1),b(2)-d(2),b(3)-d(3)]
        ca_vec=[a(1)-c(1),a(2)-c(2),a(3)-c(3)]
        v34=[gx4-gx3,gy4-gy3,gz4-gz3]
        v32=[gx2-gx3,gy2-gy3,gz2-gz3]
        
        ! 计算归一化的板坐标轴
        lbx = v34
        lbx = lbx/norm2(lbx)
        lbz = cross_product(v34,v32)
        lbz = lbz/norm2(lbz)
        lby = cross_product(lbz,lbx)

        ! 计算板坐标系的坐标值
        x3b=0.0_real64
        y3b=0.0_real64
        y4b=0.0_real64
        x4b=sqrt((gx4-gx3)*(gx4-gx3)+(gy4-gy3)*(gy4-gy3)); 
        ! 计算node 1
        g1=[gx1,gy1,gz1]
        g2=[gx3,gy3,gz3]
        call global_to_local_2d(g1,g2,lbx,lby,x1b,y1b)
        ! 计算node 2
        g1=[gx2,gy2,gz2]
        call global_to_local_2d(g1,g2,lbx,lby,x2b,y2b)
        
        ! 计算系数Ax,Bx...
        ax=x1b-x2b-x3b+x4b
        ay=y1b-y2b-y3b+y4b
        bx=x1b-x2b+x3b-x4b
        by=y1b-y2b+y3b-y4b
        cx=x1b+x2b-x3b-x4b
        cy=y1b+y2b-y3b-y4b

        ! 计算bd_vec,ca_vec的分别和lbx的家教alpha,beta的三角函数值
        call vector_angle(bd_vec,lbx,ca,sa)
        call vector_angle(ca_vec,lbx,cb,sb)
        t1(1,:)=[sb,-sa]
        t1(2,:)=[-cb,ca]

        ! 将epsilon_{xi z}填充到t2(2,12)的第一行
        t2(1,1)=0.5_real64*(1.0_real64+eta)               ! w1
        t2(1,2)=-0.25_real64*(1.0_real64+eta)*(y1b-y2b)    ! theta1_x
        t2(1,3)=0.25_real64*(1.0_real64+eta)*(x1b-x2b)    ! theta1_y
        t2(1,4)=-0.5_real64*(1.0_real64+eta)              ! w2
        t2(1,5)=-0.25_real64*(1.0_real64+eta)*(y1b-y2b)   ! theta2_x
        t2(1,6)=0.25_real64*(1.0_real64+eta)*(x1b-x2b)    ! theta2_y
        t2(1,7)=-0.25_real64*(1.0_real64-eta)             ! w3
        t2(1,8)=-0.25_real64*(1.0_real64-eta)*(y4b-y3b)   ! theta3_x
        t2(1,9)=0.25_real64*(1.0_real64-eta)*(x4b-x3b)    ! theta3_y6
        t2(1,10)=0.5_real64*(1.0_real64-eta)              ! w4
        t2(1,11)=-0.25_real64*(1.0_real64-eta)*(y4b-y3b)   ! theta4_x
        t2(1,12)=0.25_real64*(1.0_real64-eta)*(x4b-x3b)   ! theta4_y
        t2(1,:)=(0.125_real64*sqrt((cx+xi*bx)*(cx+xi*bx)+(cy+xi*by)*(cy+xi*by))/&
                det_jacobian(xi,eta,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4))*t2(1,:)
        ! ! 同理
        t2(2,1)=0.5_real64*(1.0_real64+xi)                ! w1
        t2(2,2)=-0.25_real64*(1.0_real64+xi)*(y1b-y4b)    ! theta1_x
        t2(2,3)=0.25_real64*(1.0_real64+xi)*(x1b-x4b)     ! theta1_y
        t2(2,4)=0.5_real64*(1.0_real64-xi)                ! w2
        t2(2,5)=-0.125_real64*(1.0_real64-xi)*(y2b-y3b)   ! theta2_x
        t2(2,6)=0.25_real64*(1.0_real64-xi)*(x2b-x3b)     ! theta2_y
        t2(2,7)=-0.5_real64*(1.0_real64-xi)               ! w3
        t2(2,8)=-.25_real64*(1.0_real64-xi)*(y2b-y3b)     ! theta3_x
        t2(2,9)=0.25_real64*(1.0_real64-xi)*(x2b-x3b)     ! theta3_y
        t2(2,10)=-0.5_real64*(1.0_real64+xi)              ! w4
        t2(2,11)=-0.25_real64*(1.0_real64+xi)*(y1b-y4b)   ! theta4_x
        t2(2,12)=0.25_real64*(1.0_real64+xi)*(x1b-x4b)    ! theta4_y
        t2(2,:)=(0.125_real64*sqrt((ax+eta*bx)*(ax+eta*bx)+(ay+eta*by)*(ay+eta*by))/det_jacobian(xi,eta,sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4))*t2(2,:)
        ! t1*t2得到剪切应变矩阵
        bs=matmul(t1,t2)
    end subroutine get_shear_bmatrix

    !> 膜内的应力应变矩阵
    function dm(E,nu) result(D)
        implicit none
        real(real64), intent(in) :: E,nu
        real(real64) :: D(3,3),a
        D = 0.0_real64
        a=E/(1.0_real64-nu*nu)
        D(1,1)=a
        D(2,2)=a
        D(1,2)=nu*a
        D(2,1)=nu*a
        D(3,3)=0.5_real64*(1.0_real64-nu)*a
    end function dm

    !> 板弯曲的应力应变矩阵
    function d_bend(e,t,nu) result(d)
        implicit none
        real(real64), intent(in) :: t,nu,e
        real(real64) :: d(3,3),a
        d(:,:) = 0.0_real64
        a=(e/(1-nu*nu))*((t*t*t)/12.0_real64)
        ! a=((t*t*t)/12.0)
        d(1,1)=a
        d(2,2)=a
        d(1,2)=nu*a
        d(2,1)=nu*a
        d(3,3)=0.5*(1-nu)*a
    end function d_bend

    !> 板的剪切应力-应变矩阵
    function d_shear(e,nu,t) result(d)
        implicit none
        real(real64), intent(in) :: e,nu,t
        real(real64) :: d(2,2),a
        d(:,:) = 0.0_real64
        a=(0.5_real64*e*t*k_)/(1.0_real64+nu)
        d(1,1)=a
        d(2,2)=a
    end function d_shear

    !> 膜单元平动自由度的形函数
    function tshape_func(i, xi, eta) result(retval)
        implicit none
        real(real64), intent(in) :: xi,eta
        integer, intent(in) :: i
        real(real64) :: retval
        select case(i)
            case(1)
                retval=0.25_real64*(1.0_real64+xi)*(1.0_real64+eta)
            case(2)
                retval=0.25_real64*(1.0_real64-xi)*(1.0_real64+eta)
            case(3)
                retval=0.25_real64*(1.0_real64-xi)*(1.0_real64-eta)
            case(4)
                retval=0.25_real64*(1.0_real64+xi)*(1.0_real64-eta)
        end select
    end function tshape_func
    
    !> 膜单元屏东自由度的形函数对xi的偏导数
    function tshape_func_pxi(i,xi,eta) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in), optional :: xi,eta
        real(real64) :: retval
        if (present(eta) .and. present(xi)) then
            select case(i)
                case(1)
                    retval=0.25_real64*(eta+1.0_real64)
                case(2)
                    retval=-0.25_real64*(eta+1.0_real64)
                case(3)
                    retval=0.25_real64*(eta-1.0_real64)
                case(4)
                    retval=0.25_real64*(1.0_real64-eta)
                end select
        end if
    end function tshape_func_pxi
    
    !> 膜单元屏东自由度的形函数对eta的偏导数
    function tshape_func_peta(i,xi,eta) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in), optional :: xi,eta
        real(real64) :: retval
        if (present(eta) .and. present(xi)) then
            select case(i)
                case(1)
                    retval=0.25_real64*(xi+1.0_real64)
                case(2)
                    retval=0.25_real64*(1.0_real64-xi)
                case(3)
                    retval=0.25_real64*(xi-1.0_real64)
                case(4)
                    retval=0.25_real64*(-1.0_real64-xi)
                end select
        end if
    end function tshape_func_peta

    !> 膜单元由节点附加角度产生的位移u_{\theta}的形函数
    function u_theta_shape_func(i,xi,eta,y1,y2,y3,y4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,y1,y2,y3,y4
        real(real64) :: ixi,ieta,b1,b2,b3,retval
        b1=0.0_real64
        b2=0.0_real64
        b3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        b1=0.25_real64*(y1+y4-y2-y3)
        b2=0.25_real64*(y1+y2-y3-y4)
        b3=0.25_real64*(y1+y3-y2-y4)
        retval = 0.125_real64*((1.0_real64-xi*xi)*ixi*(b1+b3*ieta)*(1.0_real64+ieta*eta)+&
                        0.125_real64*(1.0_real64-eta*eta)*ieta*(b2+b3*ixi)*(1.0_real64+ixi*xi))
    end function u_theta_shape_func
    
    !> 膜单元由节点附加角度产生的位移u_{\theta}的形函数对xi的偏导数:Ni^{u}_xi
    function u_theta_shape_func_pxi(i,xi,eta,y1,y2,y3,y4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,y1,y2,y3,y4
        real(real64) :: ixi,ieta,b1,b2,b3,retval
        b1=0.0_real64
        b2=0.0_real64
        b3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        b1=0.25_real64*(y1+y4-y2-y3)
        b2=0.25_real64*(y1+y2-y3-y4)
        b3=0.25_real64*(y1+y3-y2-y4)
        retval = 0.125_real64*(-2.0_real64*xi*ixi*(b1+b3*ieta)*(1.0_real64+ieta*eta)+&
                0.125_real64*ieta*ixi*(1.0_real64-eta*eta)*(b2+b3*ixi))
    end function u_theta_shape_func_pxi

    ! 膜单元由节点附加角度产生的位移u_{\theta}的形函数对eta的偏导数:Ni^{u}_eta
    function u_theta_shape_func_peta(i,xi,eta,y1,y2,y3,y4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,y1,y2,y3,y4
        real(real64) :: ixi,ieta,b1,b2,b3,retval
        b1=0.0_real64
        b2=0.0_real64
        b3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        b1=0.25_real64*(y1+y4-y2-y3)
        b2=0.25_real64*(y1+y2-y3-y4)
        b3=0.25_real64*(y1+y3-y2-y4)
        retval = 0.125_real64*((1.0_real64-xi*xi)*ixi*(b1+b3*ieta)*ieta+&
                        0.125_real64*(-2.0_real64*eta)*ieta*(b2+b3*ixi)*(1.0_real64+ixi*xi))
    end function u_theta_shape_func_peta

    ! 膜单元由节点附加角度产生的位移v_{\theta}的形函数
    function v_theta_shape_func(i,xi,eta,x1,x2,x3,x4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4
        real(real64) :: ixi,ieta,a1,a2,a3,retval
        
        a1=0.0_real64
        a2=0.0_real64
        a3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        a1=0.25_real64*(x1+x4-x2-x3)
        a2=0.25_real64*(x1+x2-x3-x4)
        a3=0.25_real64*(x1+x3-x2-x4)
        retval = -0.125_real64*((1.0_real64-xi*xi)*ixi*(a1+a3*ieta)*(1.0_real64+ieta*eta)+ &
                        0.125_real64*(1.0_real64-eta*eta)*ieta*(a2+a3*ixi)*(1.0_real64+ixi*xi))
    end function v_theta_shape_func

    ! 膜单元由节点附加角度产生的位移v_{\theta}的形函数对xi的偏导数:Ni^{v}_xi
    function v_theta_shape_func_pxi(i,xi,eta,x1,x2,x3,x4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4
        real(real64) :: ixi,ieta,a1,a2,a3,retval
        
        a1=0.0_real64
        a2=0.0_real64
        a3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        a1=0.25_real64*(x1+x4-x2-x3)
        a2=0.25_real64*(x1+x2-x3-x4)
        a3=0.25_real64*(x1+x3-x2-x4)
        retval = -0.125_real64*((-2.0_real64*xi)*ixi*(a1+a3*ieta)*(1.0_real64+ieta*eta)+ &
                        0.125_real64*(1.0_real64-eta*eta)*ieta*(a2+a3*ixi)*ixi)
    end function v_theta_shape_func_pxi

    ! 膜单元由节点附加角度产生的位移v_{\theta}的形函数对xi的偏导数:Ni^{v}_eta
    function v_theta_shape_func_peta(i,xi,eta,x1,x2,x3,x4) result(retval)
        implicit none
        integer, intent(in) :: i
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4
        real(real64) :: ixi,ieta,a1,a2,a3,retval
        
        a1=0.0_real64
        a2=0.0_real64
        a3=0.0_real64
        ixi =node_natural_coords(i,1)
        ieta=node_natural_coords(i,2)
        
        a1=0.25_real64*(x1+x4-x2-x3)
        a2=0.25_real64*(x1+x2-x3-x4)
        a3=0.25_real64*(x1+x3-x2-x4)

        retval = -0.125_real64*((1.0_real64-xi*xi)*ixi*(a1+a3*ieta)*ieta+ &
                        0.125_real64*(-2.0_real64*eta)*ieta*(a2+a3*ixi)*(1.0_real64+ixi*xi))
    end function v_theta_shape_func_peta

    ! 雅可比矩阵J eq2.28 p20
    function f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4) result(jacb)
        implicit none
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4,y1,y2,y3,y4
        real(real64) :: jacb(2,2)
        real(real64) :: x_xi,x_eta,y_xi,y_eta
        x_xi    = 0.25_real64 * ((1.0_real64+eta)*x1+(-1.0_real64-eta)*x2+(eta-1.0_real64)*x3+(1.0_real64-eta)*x4)
        x_eta   = 0.25_real64 * ((1.0_real64+xi)*x1 +(1.0_real64-xi)*x2+(xi-1.0_real64)*x3+(-1.0_real64-xi)*x4)
        y_xi    = 0.25_real64 * ((1.0_real64+eta)*y1+(-1.0_real64-eta)*y2+(eta-1.0_real64)*y3+(1.0_real64-eta)*y4)
        y_eta   = 0.25_real64 * ((1.0_real64+xi)*y1+(1.0_real64-xi)*y2+(xi-1.0_real64)*y3+(-1.0_real64-xi)*y4)
        jacb(1,1) = x_xi
        jacb(2,1) = x_eta
        jacb(1,2) = y_xi
        jacb(2,2) = y_eta
    end function f_jacobian

    ! 雅可比矩阵J的行列式
    function det_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4) result(det)
        implicit none
        real(real64), intent(in) :: xi,eta,x1,x2,x3,x4,y1,y2,y3,y4
        real(real64) :: det
        real(real64) :: jacb(2,2)
        jacb=f_jacobian(xi,eta,x1,x2,x3,x4,y1,y2,y3,y4)
        det=jacb(1,1)*jacb(2,2)-jacb(1,2)*jacb(2,1)
    end function det_jacobian
    
    subroutine vector_angle(a,b,cos_theta,sin_theta)
        implicit none
        real(real64) ,dimension(3),intent(in) :: a, b
        real(real64),intent(out) :: cos_theta
        real(real64),optional,intent(out) :: sin_theta
        real(real64):: norm_a, norm_b, dot_ab, norm_cross
        ! 计算点积 a·b
        dot_ab = dot_product(a, b)
        ! 计算向量模长 ||a|| 和 ||b||
        norm_a = norm2(a)
        norm_b = norm2(b)
        ! 计算余弦值
        cos_theta = dot_ab / (norm_a * norm_b)
        ! 计算叉积模长 ||a×b||（3维向量专用）
        norm_cross = norm2([a(2)*b(3) - a(3)*b(2), &
                            a(3)*b(1) - a(1)*b(3), &
                            a(1)*b(2) - a(2)*b(1)])
        
        if (present(sin_theta)) then
            ! 计算正弦值
            sin_theta = norm_cross / (norm_a * norm_b)
        end if
    end subroutine vector_angle

    function cross_product(a,b) result(c)
        ! ! deepseek
        implicit none
        real(real64), dimension(3), intent(in) :: a, b
        real(real64), dimension(3) :: c
        c =[a(2)*b(3) - a(3)*b(2), &
            a(3)*b(1) - a(1)*b(3), &
            a(1)*b(2) - a(2)*b(1)]
    end function cross_product

    subroutine global_to_local_2d(A_global, C_global, x_axis, y_axis, ax,ay)
        implicit none
        real(real64), intent(in) :: A_global(3)      ! 点 A 的全局坐标 (A_x, A_y, A_z)
        real(real64), intent(in) :: C_global(3)      ! 局部坐标系原点 C 的全局坐标 (C_x, C_y, C_z)
        real(real64), intent(in) :: x_axis(3)        ! 局部坐标系 x 轴基向量（已归一化）
        real(real64), intent(in) :: y_axis(3)        ! 局部坐标系 y 轴基向量（已归一化）
        real(real64), intent(out) ::ax,ay           ! 点 A 在局部坐标系中的坐标 (A_x, A_y)
        logical:: is_in_plane  ! 是否在平面内（否则投影）
        real(real64) :: r(3), n(3), r_proj(3), r_dot_n
        ! 计算相对位移向量 r = A - C
        r = A_global - C_global
        ! 计算法向量 n = x × y
        n(1) = x_axis(2) * y_axis(3) - x_axis(3) * y_axis(2)
        n(2) = x_axis(3) * y_axis(1) - x_axis(1) * y_axis(3)
        n(3) = x_axis(1) * y_axis(2) - x_axis(2) * y_axis(1)
        ! 检查 A 是否在 x-y 平面内（r · n ≈ 0）
        r_dot_n = dot_product(r, n)
        is_in_plane = (abs(r_dot_n) < 1e-6)  ! 浮点误差容忍度
        ! 计算投影向量 r_proj = r - (r · n) n
        if (is_in_plane) then
            r_proj = r
        else
            r_proj = r - r_dot_n * n
        end if
        ! 计算局部坐标 A_X = r_proj · x, A_Y = r_proj · y
        ax = dot_product(r_proj, x_axis)
        ay = dot_product(r_proj, y_axis)
    end subroutine global_to_local_2d

    subroutine global_to_local(a_global, origin_B, x_axis_B, y_axis_B, z_axis_B, a_local)
        !> @brief: 计算空间点a坐标值转换到在局部坐标的坐标值
        !> 
        !> @param[in] a_global: 全局坐标系中的点A坐标 (ax,ay,az)
        !> @param[in] origin_B: 局部坐标系B的原点全局坐标
        !> @param[in] x_axis_B: 局部坐标系B的x轴单位向量 (已归一化)
        !> @param[in] y_axis_B: 局部坐标系B的y轴单位向量 (已归一化)
        !> @param[in] z_axis_B: 局部坐标系B的z轴单位向量 (已归一化)
        !> @param[out] a_local: 输出: A点在局部坐标系B中的坐标 (AX,AY,AZ)
        !> 
        !> @note: 局部坐标系B的坐标轴定义为:
        !>         - x轴: 局部坐标系B的x轴
        !>         - y轴: 局部坐标系B的y轴
        !>         - z轴: 局部坐标系B的z轴
        !>         - 原点: 局部坐标系B的原点
        !>         - 基向量: 局部坐标系B的x轴、y轴、z轴
        implicit none
        real(real64), intent(in) :: a_global(3)       ! 全局坐标系中的点A坐标 (ax,ay,az)
        real(real64), intent(in) :: origin_B(3)       ! 局部坐标系B的原点全局坐标
        real(real64), intent(in) :: x_axis_B(3)      ! 局部坐标系B的x轴单位向量 (已归一化)
        real(real64), intent(in) :: y_axis_B(3)      ! 局部坐标系B的y轴单位向量 (已归一化)
        real(real64), intent(in) :: z_axis_B(3)      ! 局部坐标系B的z轴单位向量 (已归一化)
        real(real64), intent(out) :: a_local(3)       ! 输出: A点在局部坐标系B中的坐标 (AX,AY,AZ)
        real(real64) :: relative_pos(3)  ! 相对位置向量
        ! 1. 计算点A相对于局部坐标系原点的位移
        relative_pos = a_global - origin_B
        ! 2. 通过点积计算在各局部坐标轴上的投影分量
        a_local(1) = dot_product(relative_pos, x_axis_B)  ! AX分量
        a_local(2) = dot_product(relative_pos, y_axis_B)  ! AY分量
        a_local(3) = dot_product(relative_pos, z_axis_B)  ! AZ分量
    end subroutine global_to_local

    !>  输入四节点壳单元的节点全局坐标值,计算
    subroutine get_shell_sys_coords(gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,&
                                    sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4)
        real(real64), intent(in) :: gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4
        real(real64), intent(out) :: sx1,sx2,sx3,sx4,sy1,sy2,sy3,sy4
        real(real64):: global_p(3),local_p(3),shell_sys_origin(3),&
                shell_sys_x(3),shell_sys_y(3),shell_sys_z(3)
        ! 计算节点在壳单元局部坐标系下的坐标值
        global_p(:)=0.0_real64
        shell_sys_origin=0.25_real64*[gx1+gx2+gx3+gx4, gy1+gy2+gy3+gy4, gz1+gz2+gz3+gz4]
        shell_sys_x=0.5_real64*[gx1+gx4, gy1+gy4, gz1+gz4]-shell_sys_origin
        shell_sys_y=0.5_real64*[gx1+gx2, gy1+gy2, gz1+gz4]-shell_sys_origin
        shell_sys_z=cross_product(shell_sys_x,shell_sys_y)
        shell_sys_x=shell_sys_x/norm2(shell_sys_x)
        shell_sys_y=shell_sys_y/norm2(shell_sys_y)
        shell_sys_z=shell_sys_z/norm2(shell_sys_z)
        
        global_p=[gx1,gy1,gz1]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx1=local_p(1)
        sy1=local_p(2)
        
        global_p=[gx2,gy2,gz2]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx2=local_p(1)
        sy2=local_p(2) 
        
        global_p=[gx3,gy3,gz3]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx3=local_p(1)
        sy3=local_p(2)
        
        global_p=[gx4,gy4,gz4]
        call global_to_local(global_p,shell_sys_origin,shell_sys_x,shell_sys_y,shell_sys_z,local_p)
        sx4=local_p(1)
        sy4=local_p(2) 
    end subroutine get_shell_sys_coords

    !> 计算位移坐标转换矩阵,U_local=T*U_global
    subroutine get_tmatrix(gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4,tmatrix)
        use matrix_utils_mod , only : invert_matrix_d
        implicit none
        real(real64), intent(in) :: gx1,gx2,gx3,gx4,gy1,gy2,gy3,gy4,gz1,gz2,gz3,gz4
        real(real64), intent(out) :: tmatrix(24,24)
        real(real64):: global_sys(3,3),shell_sys_origin(3),&
                shell_sys(3,3),temp(3,3),temp_inv(3,3)
        integer:: i,j
        global_sys(:,1)=[1.0_real64 , 0.0_real64 , 0.0_real64]
        global_sys(:,2)=[0.0_real64 , 1.0_real64 , 0.0_real64]
        global_sys(:,3)=[0.0_real64 , 0.0_real64 , 1.0_real64]

        shell_sys_origin=0.25_real64*[gx1+gx2+gx3+gx4, gy1+gy2+gy3+gy4, gz1+gz2+gz3+gz4]
        shell_sys(:,1)=0.5_real64*[gx1+gx4, gy1+gy4, gz1+gz4]-shell_sys_origin
        shell_sys(:,2)=0.5_real64*[gx1+gx2, gy1+gy2, gz1+gz4]-shell_sys_origin
        shell_sys(:,3)=cross_product(shell_sys(:,1),shell_sys(:,2))
        shell_sys(:,1)=shell_sys(:,1)/norm2(shell_sys(:,1))
        shell_sys(:,2)=shell_sys(:,2)/norm2(shell_sys(:,2))
        shell_sys(:,3)=shell_sys(:,3)/norm2(shell_sys(:,3))
        
        do j=1,3
            do i=1,3
                temp(i,j)=dot_product(shell_sys(:,j),global_sys(:,i))
            end do
        end do
        
        call invert_matrix_d(temp,temp_inv,j)
        if (j /= 0) then
            print *, "Error happened at subroutine : get_tmatrix"
        end if

        tmatrix=0.0_real64
        do i=1,8
            tmatrix([3*i-2,3*i-1,3*i],[3*i-2,3*i-1,3*i])=temp_inv
        end do
    end subroutine get_tmatrix
end module shell_mod